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Abstract 



Spatial non-locality of space-fractional viscoelastic equations of motion 
is studied. Relaxation effects are accounted for by replacing second-order 
time derivatives by lower-order fractional derivatives and their generaliza- 
tions. It is shown that space-fractional equations of motion of an order 
strictly less than 2 allow for a new kind anisotropy, associated with an- 
gular dependence of non-local interactions between stress and strain at 
, different material points. Constitutive equations of such viscoelastic me- 

. dia are determined. Explicit fundamental solutions of the Cauchy problem 

QQ ' are constructed for some cases isotropic and anisotropic non-locality. 

I Keywords, viscoelasticity, non-local, fractional, anisotropic, bio-tissue, anoma- 

' ' lous diffusion. MSG Class: 74D05, 74A20, 74J10. 

Notation. 

(/** gm ■.^J^fis)g{t-s)ds 
if *x 9)ix) := /jj3 f{x - y) g{y) Ay 
Fourier transform: /(k) — J^^ Q-ik x j^^^ 
Laplace transform: f{p) — e^^* f{t) dt. 

rS 

H . 

- - ■ 1 Introduction 

Partial differential equations with fractional derivatives play an important role 
in modeling anomalous diffusion and wave propagation in media with complex 
structure. Equations invariant with respect to space and time scaling are of 
particular interest. While space-time fractional equ ations have be en widely 



accepted and studied for modeling anomalous diffusion iHanygal [2002|, no anal- 
ogous equations have been obtained for wave propagation in continuous media. 

Viscoelastic models considered here have two new aspects. In the first place 
the stress-strain constitutive equation is non-local - the stress at a material 
point depends on the strain at neighboring material points - and as a result 
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the operator in the equation of motion acting on the spatial variables is non- 
local. In addition, the time derivative can be replaced by a non-local operator 
to allow for time-delayed response of stress to strain. The latter addition is 
however familiar in contemporary viscoelasticity. 

The other aspect of spatial non-locality is the appearance of a new kind of 
anisotropy. In linear viscoelasticity elastic anisotropy manifests itself in connec- 
tion with the tensorial character of the constitutive equation linking stress to 
the strain tensor. Elastic anisotropy is associated with the symmetry properties 
of the stiffness tensor with respect to the linear transformations of the coordi- 
nates. Elastic anisotropy can result from crystallographic structure, fine layered 
structure, parallel cracks etc. It accounts for angular dependence of elastic and 
viscoelastic moduli, that the dependence of stress on the direction of stretching 
or shear. 

In the theory presented below there is room for a different kind of anisotropy 
associated with non-locality of the stress-strain constitutive equation. In the 
equations of motion the operator acting on the spatial variables is a pseudo- 
differential operator of an order < 2. If the order of the spatial operator is 
strictly lower than 2, then the spatial operator operator is non-local. This 
brings about a new kind of anisotropy. The new anisotropy is not related to 
the tensorial properties of the stiffness coefficients but is instead defined by 
the micro-local structure of the spatial operator. Anisotropy associated with 
non-locality accounts for the angular dependence of sensitivity of stress to the 
strain in the adjacent material. Such a dependence can be expected to be much 
stronger along a ligament than transversally to it. 

In the Fourier-transformed spatial operator the anisotropy associated with 
non-locality is expressed in terms of anisotropic wave number dependence. We 
shall only consider a special class of such operators whose symbols are expressed 
in terms of a probability distribution of over all the directions on the unit sphere. 
A distribution over directions on a unit sphere can be expressed as an infinite 
series of Legendre polynomials. This is far more general than the concept of 
anisotropy considered in elasticity. 

In linear elasticity all the anisotropy classes are derived from the symmetries 
of the stiffness coefficients. The anisotropic properties of such media can be 
expressed in terms of tensors of finite rank. In local linear viscoelasticity stiffness 
coefficients are frequency dependent. Such understanding of anisotropy has 
been motivated by crystallographic, layered or densely cracked structures of 
many materials. In modelling of high-resolution MRI scans of anisotropic bio - 
tissues requires diffusion tensors of very high orders [Ozarslan and Mareci 2003| . 
This is due to a high complexity of cell and tissue structure. Anisotropy of 
viscoelastic wave motion in a bio-tissue can be expected to reveal the same level 
of complexity. Anisotropic diffusion is possible only for anomalous diffusion of 
the space-fractional type. Anomalous diffusion has a counterpart in non-local 
viscoelasticity of the kind considered in this paper. 

The anisotropy of non-locality can coexist with the anisotropy associated 
with the symmetries of the stiffness coefficients because the stiffness coefficients 
depend on the wave number directions. This kind of anisotropy is however also 



2 



possible in scalar models of viscoelasticity. It is interesting that for this kind 
of anisotropy explicit solutions of quite general anisotropic equations can be 
constructed, which is not the case for anisotropy associated with the symmetries 
of the stiffness coefficients. 

In the case of anisotropy associated with the symmetries of the stiffness coef- 
ficients the simplest geometrical element is a plane of symmetry (in transversal 
symmetry). Planes of symmetry can be associated with micro-layered struc- 
tures. Anisotropy associated with non-locality allows for different kinds of mi- 
crostructure, such as curvilinear ligaments and channels in bio-tissues. Non- 
locality is also very likely in bio-tissues due to their complex structure and 
composition. 

In this paper we focus on spatial pdo's whose symbols are homogeneous 
functions of the wavenumber. In the absence of ani s otrop y such operators reduce 
to fractional-order Riesz derivatives Samko et al. 1993l | or, in other words, to 
fractional-order Laplacians. 

The usual Laplace-Fourier representations are inconvenient for the construc- 
tion and analysis of the solutions of equations which are of fractional order 
with respect to both spatial variables and time variable. This leads to some 
difficulties in the ana lysis of the attenuation and dispersion. As discovered in 
Mainardi et al. I (200l| . the solutions of such equations can be expressed in terms 



of Mellin convolutions. The solution is a superposition of copies of a fixed wave- 
form subject to a varying scaling. The above-mentioned representations of the 
solutions are fairly explicit and involve only one non-elementary function - the 
Wright function. The Wright function corrects for the difference between the 
orders of the time derivative and the spatial operator. The effect of the non- 
local spatial operator is represented by an elementary function provided the the 
order of the time derivative matches the order of the Laplacian. 



2 Space-time fractional wave equations 

We shall derive an explicit formula for the Green's function of an anisotropic 
space-time fractional linear viscoelastic equation, defined as the solution of the 
Cauchy problem 

pB^^u^Qu, u{0,x)^S{x), Du(0,x)=0 xeR^ (1) 

where denotes the Caputo derivative, 1 < /? < 2 and Q is a pseudo- 
differential operator defined by its symbol 

Q(k) = -^|k.yrMdy) (2) 

and (3 < a < 2. The integral in eq. ([2]) extends over a unit sphere S, defined 
by the equation |y| = 1 and the measure fi on S is non- negative and has finite 
mass. 

Note that Q(k) — |k|" Q(k). Hence the order a of the operator is indepen- 
dent of the anisotropic properties of Q, represented by the measure /u. 
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Choosing /i to be a homogeneous distribution over the sphere 

a + 1 

^(dy) = M — — sin(6i) d9 dip 
where 0, ip are polar coordinates on S we get 

Q(k) = -M|kr (3) 
and thus Q essentiahy reduces to a fractional-order Laplacian 

Q = -M{-S/Y' (4) 

For a = 2 we have |k • yp = k • yy^ y and thus Q(k) = — k • M k, where M 
is a positive definite matrix 



M:= / yy^M(dy) 
Js 



Consequently Q — V • MV is an anisotropic generalization of the Lapla- 
cian. A spatial operator of second order is ellipsoidal allows only for ellipsoidal 
anisotropy. 

A different anisotropic generalization of the Laplacian is obtained by choos- 
ing a measure fi with three mutually orthogonal support points yi , y2 , ya on 
the sphere S. In this case Q(k) = —M |k|" k • Ak, where k |k|^^ k. 

Let P denote the fractional integral operator 

f* (t ~ sV^^ 

^ p(^) /G^)d., ^>0 (5) 

for 7 > 0. For positive integer /3 the Caputo derivative is an ordinary derivative. 
For positive non-integer /? Caputo derivative D'^ of order /3 is defined by the 
formula 

D'3/ = I"-'3d"/ (6) 



where n is an integer such that n — 1 < /3 < n iPodlubnvi [1998 1. 



3 Comparison with scalar viscoelasticity 

3.1 Construction of fractional-order viscoelastic equations 
of motion. 

As it stands, equation ([1]) does not have the form of a viscoelastic equation 
of motion. We shall prove that it is equivalent to a generic momentum bal- 
ance equation p m = div D cr, where the stress cr is given by the constitutive 
equation 

= *t CDe, e = (Vm + (Vu)"^) /2 (7) 
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g is a completely monotone function (appendix [X]) and C is a non-local integral 
operator of order a — 2 acting on the space variables x. 

We shall begin with transforming equation ([l} into a second order equation 
with respect to time. 

Equation ([T]) implies that u{t, x) — u(0, x) — t u{0, x) = B'^ u = p"^ l^^ Q u, 
hence pD^u = D^l^Qu. The right-hand side can be expressed in the form 
DI'5-iQu = DI'5-iQDu+ [t^-Vr(/3)] Qu{0,x). 

pVp- u^Vp-\^ Qu (8) 

The right-hand side of equation ^ can be transformed as follows: 

D2 f Tf>-'/T{P){Qu){t-T)dT^ 

Jo 

Dt'^-Vr(/3)(Qu)(0)+D / T^-'/r{f3){QBu){t-T)dT = 

Jo 

Dtl-'/T{fi)H {QDu){t) = 

f\t-Tf-yr{/3-l){QDu){T)dT (9) 

where we have taken advantage of the inequality P > I and the commutation of 
Q and D and set (Du){0,x) = 0. We thus have transformed equation ^ into 
the following form 

pB^ u = g*t Qu (10) 

where the relaxation modulus g{t) = t^~'^ /r{P — 1). The relaxation modulus g 
is locally integrable completely monotone (LICM, see Appendix \K\ 

In the next step we shall recast the equation in the form of a momentum 
conservation equation by setting Qu = divK(it). To this effect we note that Q, 
given by equation is a homogeneous and difFerentiable function of degree a: 

Q(Ak) - A"Q(k) 

Hence Q(k) = ik • K(k), where K = Vk Q(k) / (ia). We now define the operator 
K by its symbol K. The stress (momentum flux) can now be identified as 
cr = (7*f D K(u) and equation (jlOp assumes the form of a momentum conservation 
equation: 

pD2w = div(j (11) 

We now verify whether the stress is a functional of the strain rate e. To this 
effect we note that (3(k) = -k-C(k) k, where C(k) := -Vk VkQ(k)/[a (a-1)]. 
Hence K(k) = i C(k) k/(a - 1), K(u) = C(-iV)w and a = .g *t D C(-iV)M = 
g*t Q(De), where Q is a non-local operator acting on spatial variables. Conse- 
quently stress is given as a non-local linear functional of the strain history. C 
is the extension of the stiffness tensor to the spatially non-local viscoelasticity. 
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The above procedure can be readily extended to vectorial equations of vis- 
coelasticity. In this case Q is a tensor- valued operator 



(Qu)fe - Qki ui (12) 
where u denotes the displacement vector, and 

Cfei = .9 *t Q/c;mn(DUm,n) 

We now require that (1) the stress effectively depends only on strain e = 
(Vu + (Vu)^) /2, (2) symmetry of stress a = follows from directly from 
the constitutive equations. The two requirements will be satisfied if Ckimn = 
Ckinm = Cikmn- The abovc symmetries are satisfied if there is a homogeneous 
generating function V^(k) of degree a + 2, with continuous derivatives up to 4-th 
order, such that 

For the operator defined by ^ the generating function assumes the following 
form 

^^^^ = ^r^TTT^ / 1"" • ^'^'^y^ ^i^) 

(Q + l)(a + 2) Js 

For simplicity we have assumed here that the kernel g is scalar. Otherwise 
we would have to take into account commutation of the tensor- valued functions 
or else to consider convolution with a tensor-valued kernel dependent on both 
time and space variables. Viscoelastic equations with tensor-valued relaxation 



modu les g and local spatial operators are considered in lHanyga and Seredynska 
120071. 



3.2 Energy conservation and its implications. 

Existence of an energy conservation with a non-negative energy functional plays 
an important role in proofs of well-posedncss of the equations. We shall examine 
some assumptions on the operator Q that ensure existence of a non-negative 
energy satisfying a conservation equation. 

The following energy balance holds for solutions of equation ^ (3 — 2 with 
Q — div o C o grad : 



d 
dt 



/ ^pii^dx + l [ {\7u f C (Vw) dx 



= 



(15) 



where ii Dm. The pseudo-differential operator C can be expressed as a 
convolution with a spatial distribution H: Cf = H *x S ■ The Fourier transform 
H of the kernel H is the symbol of the operator C. The second term in the 
square brackets on the left-hand side represents the stored energy U{t) at time 
t. 
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We now impose the condition that the stored energy is non- negative. Apply- 
ing the Fourier transformation to U we have on account of Plancherel's theorem 

U{0) = 1 [ (ikw(0,k))^i7(k)(ikw(0,k))dkEE~i / u(0, k)^ Q(k) ti(0, k) dk 

(16) 

In view of the arbitrariness of the function u(0, x) we shall impose the condition 
that Q(k) is Hermitian positive semidefinite 

g(k) < (17) 

i.e. w'f (5(k) w for all complex vectors w e C'^. 

In particular, for scalar equations with Q = —M (— V^)"^^, 1 < a < 2, we 
have g(k) = -M (k^)"/^ < 0, and 

because 2 — a < 1. 

Let us now turn to 1 < (3 <2. The equation p T)^ u = Qu can be transformed 
to the more familiar form 

If the stress is defined by the equation 

a = g*tH *x VDu (19) 

then equation (fTO|) assumes the familiar form ([TT]) . Taking the scalar product 
of both sides of the last equation with the vector Du we arrive at the formula 

A^(Dw)2 = -trace (ffVDu) (20) 

We shall now construct the stored energy functional U in such a way that 
dU/dt = trace (cj VDw). 

The function g{t) is completely monotone and lo cally integrable (LICM) . 
LICM functions are causally positive definite (CPD) iGripenberg et al. 1990| |. 
Wc shall therefore carry out our analysis for a general causally positive definite 
function G. A CPD function g can be expressed as the Fourier transform of of 
a positive Radon measure m 

/oo 
e'^*m(ds) (21) 
-oo 

where 

m(ds) = -3ig{s) (22) 
In particular, for g(t) = 4~ ^/r(^-l) we have m(ds) = (I/tt) sin(^7t/2) |s|1-^ ds 



Gel'fand and Shilovl [1964 1. 
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Let Tp := H *^ VDu, 



yit,x;s) 



We now make the following assumption 

Assumption 1 For each k the matrix H(\i) is positive definite and symmetric. 

The function H an even function of k because H is real-valued. Define the 
kernel 

1 



e 



— ik x £ri,— 1 



i/k-^ dk 



(23) 



It is easy to check that Hi{x) *x H{x) — 5{x). We now define the stored energy 
functional 



U 



I 



y(t, x; s)^ Hi {x ~ y) y{t, y; s) m(ds) dy dy (24) 



Since Dy = ip — isy and is even and symmetric 
dU f°° 1 



dt 



— OO 

OO { poo 



— oo y.'J — OO 



{ip{t, x)'^ H_i{x - y) y(t, x; s) + yt H_i{x - y) ■)p{t, x)) dy dx m{ds) 



y{t,x;s)m{ds) 



H-iix -y)i;(t,y)dy\- dx 

a J 

tra.ce{(j{t,x)VDu{t,x))dx (25) 



as required. 

We have thus proved energy conservation 



d_ 

di 



■dx + U 







(26) 



with the stored energy U defined by equation ([24|) . Assumption [T] played a key 
role in the construction of the stored energy potential. 

Assumption [T] also ensures that U >0. Indeed, by Plancherel's theorem, 

C/= / y(t,k;s)t7?(k)-iy(t,k;s)m(ds)>0 

4 Solution of the Cauchy problem for 1 < < 

a<2. 

4.1 Formulation of the Cauchy problem. 

Consider an abstract time-fractional equation 

D^u^Au (27) 
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where A is an operator and D'' is the Caputo fractional derivative. For fS < 
1 the Caputo derivative — D and, using the semigroup property of 
the fractional integral operators P, equation (P7)) is equivalent to the integral 
equation 

IDu^l^ Au 



or, more explicitly. 



u{t) = u{0)+l'^ A', 



We thus expect that equation (l?f)) with the initial condition u(0) — uq has a 
unique solution. 

If 1 < /3 < 2, then D-^ = P^l^ and equation (gT]) is equivalent to D^u = 
Au, or 

u{t) = u{0) +tu'{0) +11^ Au 

where u' = Du. In this case we expect that equation ([?7l) with the initial 
conditions 



m(0) = Mo, u'{0) — iiQ 

has a unique solution. 

Applying the Laplace transformation we have 



uo p' 



+ f[p) 



for < ^ < 1 and 



for 1 < /3 < 2. 



u{p) 



uqp^ ^ + iiQp^ ^ 
+ fip) 



(28) 



(29) 



(30) 



4.2 Solution of the Cauchy problem. 

The Cauchy problem for the scalar equation ([T]), ^ in 3 dimensions has an 
explicit solution in terms of an integral over a unit sphere with a n integrand in- 
volvin g Wright functions Mainardi 2011 1 , iMainardi et al.l 2010| , Gorenflo et al 
1999| . T he Wright fun c tions are numericall y computable by an integral repre- 



sentation lL"uchkol|2008l| . lLuchko et al.l|2010l |. An alternative method for numer- 



ical integration of would involve pseudo-spectral methods combi ned with a smart 
repres e ntation of the Caputo f r actional derivativ e the method of lYuan and Agrawal 



epres e ntation or tne >^aputo i r actionai derivativ e tne metnoa oi l Yuan ana Agrav 
1998i . ILu and HanTeal |2005| . iDiethelml |2010t . On the other hand, for some 



values of the index the Wright function can be expressed in terms of the expo- 
nential and the Airy function. 

Applying to equation ([T]) the Laplace transform f{t) — >• £[f]{p) = f{p) = 
-C° -^^^^ Fourier trans form ^(a:) ^ Trr^m^-i = — f o-ik a; , 

and the identitv IPodlubnvl |l998i 



J-[<?](k)^5(k)=43 e-*-ff(x)dx 



£ [D^/] ip) = / f{p) - p^-i /(O) - /'(O) 



(31) 
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valid for 1 < /3 < 2, wc obtain the equation 

u{p, k) T [c[u]] {p, k) = ^^oik) + Mk) f I (32) 

pf^ + Q(k) pf^ + Q(k) 

The solution of the Cauchy problem ([1]) , ([25)1 will be obtained by inversion 
of the Laplace and Fourier transformations, 

"(i.^) = :5^ / dpeP*-^ / dke"'-[/(p,k) 

Hence 

u{t,x)^ Gp^qit^x - y)uo{y)dy + Hp^Q{t, x - y) uo{y) dy (33) 



where 

Gp.Qit.x)^^ j dpe^'-^J e^-^-^^dgk (34) 

2711 7g (271)-^ + Q(k) 

Hp,Q{t.x)^^ f dpe^*-l^ / e''-- /^' dsk (35) 
2m (27t)'^ Jj^s pP + g(k) 

and dk = dk x sm{-&) di} x dtp. 

The functions G/3^q and ff^,Q will be called the first fundamental solution 
and the second fundamental solution of the Cauchy problem for ([T]), respectively. 

Exchanging the order of the inverse Laplace and Fourier transformations we 
have 

Gp,Qit,x)^-^ [ e"'--i?;3(Q(k)i)d3k (36) 

H^Mt^ x) = j:^ I e*- Ep,^{Q{k) t) dgk (37) 

If Q = — a|k|", a > 0, then [/(p, k) is a function F(p, fc) of fc := |k| and 
the integration over the angular coordinates in (p4|) and ([35l) can be carried out 
explicitly, yielding G^'^q := G^^g, where we use a new parameter 7 := j3/a for 
reasons that will soon become clear: 

G^'Ut^x) = — / dpeP*-— ^ / fce''='-F(p,|fc|)dfc = 
2711^5 (27t)^iry_oo 

2^ h^^''^r8-rYn ^e^^., 1^1) dfc ^ -^^G^ (,0U,, (38) 

where G^^a denotes the solution of the same problem in one-dimensional space. 
Similarly 

= (39) 

Let / be an even differentiable function on the real line. The mapping 
/(r) — > F{r) := —f'{r)/{2nr) has two noteworthy properties: 
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(i) if / vanishes at infinity, then 

f{\x\)dx^ I /(r)dr; 



(ii) if / is unimodal with a mode at (i.e. if f'{r) < for r > 0), then F > 0. 

Consequently, if / is a probabihty density on the real line with a maximum or 
a singularity at 0, then F(|a;|) is a probability density in M'^. This observation 
is crucial for distinguishing between diffusion and wave propagation in 1 and 3 
dimensions. 

The scaling transformations allow considerable simplification of the problem 
by reducing the number of unknowns. Changing the integration variables p — >■ 
s = pt, k ^ K = k we get the scaling relation 



GifX{t,x)^t-^''Gifl{l,x/t'>) (40) 
H\i{t,x)^t'-'''H\'^l{l,x/t'') (41) 

A simple integ ral representation of th e G'^^L(1, x) is constructed in the paper 
of Mainardi et al. iMainardi et al.l 200 1| : 

G« (1, x) = M,(0 X^{x/0 f (42) 

where 

f i |yr~' sin(a7t/2) for < a < 2 

\(<5(y-l) + % + l))/2 fora = 2 

We derive it in a slightly different way in Appendix [C] and Appendix [E] It is 
easy to see that Xa[y) > for — oo < y < oo and 



/oo 
Xc,{y)dy = l 
-oo 



hence Xg is a probab ility distribution. Xa is the function denoted by in 
iMainardi et al.l 2001 1. The function Mj is a special case of th Wright function 
(Appendix |B]): 

M^{z) = W^,i_^(-z), < 7 < 1 (44) 



Consequently 



G(3)^(l,x) = [/(^■")(r) r M,{S,)Xl^\r/0 ^ (45) 

where r — \x\ and 7 = /3/a, while 

, 19 w"^^ sin(a7T/2) , , 

X^^\y) -.^ -—r- — - f , ,L , , 2/>0, 0<a<2 (46 

" 2t? ydyl + 2y" cos(a7t/2) + y^'^ ^ ^ ' 
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Figure 1: Scaled plots of the function Za{y) ■— " Xa (y)- (1) solid line: 
30 Zi, (2) dashed: 2OZ1.5, (3) dot-dashed: Z1.9. 

It is easy to see that Xa \y) — y"^^ Za{y), where Za is regular at 0. 

In Fig. [1] scaled plots of Ya are shown. The function Xi is the Cauchy 

(3) 

probability density and X^ is non-negative. For a < 1 the function Xa is 

(3) 

unimodular with a singularity at y = and therefore Xa is non-negative. For 
a > 1 the maximum of Xa is shifted to the right of j/ = and therefore the 
function Xa^ changes sign. 

The parameter ^ scales the spatial coordinates of the component waves Ya . 
The function M^(^) represents the weights of various scalings of the solution. 
In Fig. [2] the function M^(z) is plotted for 7 = 1/2,1/3,2,3. The shift of the 
maximum of the Wright function for 7 > 1/2 is noteworthy. 

The functions a ^ = Ij 3 and a — 1.5 and 1.9 are shown in Fig. [3] 

and Fig. |4l The values of /3 are 1, and 1.2666. 

5 Fundamental solutions for anisotropic non-locality. 

5.1 Neutral case 1 < (3 = a < 2. 

We now assume that Q(k) is given by equation ([5]). The inverse Laplace trans- 
form of G(p,k) and ^(p,k) is Ep (Q{k)t^^ and Efi^2 (Q{k)t^y respectively 

or, equivalently, Ej3 (-/c" F(k) tf^^ and Ei3'2 (-/fc" F{k) ff^^ , respectively, where 

k := |k|, k := k'^ k and F(k) = ~Q(k). The Mittag-LefHer function Ex^^{z) is 
defined in the appendix and Ex i?A,i- 
We shall begin with the case (3 = a > 1. 
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Figure 2: The function M^{z): (1) solid line: 7 = 1/2, (2) dashed line: 7 = 1/3; 
(3) dot-dashed hne: 7 = 2/3 . 



U 




Figure 3: The function G(2/3 a)/^ ~ -'^•'^ (sohd line) and a = 1.5 (dashed 

hne). 
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u 




Y 



Figure 4: The function G)J^^ ^^/2 for a = 1.9 (solid line) and a = 1.5 (dashed 
line). 

The first fundamental solution can be obtained by inverting the Laplace and 
Fourier transform: 

G{t, x) = ^ £ edk j^A2\^E^ i^^k^ F(k) r) c'"-^ = 

-V^^^ £d2kF(k)-i/«^°°£;^(_^«) e'«[(''-*)'-™^'°']d«=:G(^)(t,x) 

where x := x. The inverse Fourier transforms of Ea (—n") and £'q,2 (~'*") 
are Xa{y) and ^^(j/), respectively (Appendix [E| , hence 

Gf)(t,.) = -_1_ / -^^^^xJi.-x^:—-\ (47) 
1'"^'^ (271)3^5 F(k)i/" <F(k)i/"y ^ ^ 

/ ^^^^rJi^-x 1 \ (48) 

1'"^'^ (27t)3ii5 F(k)i/" " \^ iF(k)i/"y ^ ' 

The subscript "1" refers to the value of 7 = P/a. 

Choosing the spherical coordinates (r, i?, 93) in the x-space in such a way that 
kk corresponds to d — Q and k • x = cos(?9), and substituting 

a rather complicated but fairly explicit expression is obtained. 
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We shall denote the functions G{t, x) and H{t, x) obtained above by ^{t, x) 
and Hf^^^{t,x). 

5.2 Anisotropic non-locality, 1 < /3 < a < 2. 

For the case /? < a we shall use the identity 

/>oo 

Ep{-v)^ M^iOEo^i-vn^i (49) 

JO 

where 7 = /3/a, see Appendix [Cl Hence, setting t — 1 for simplicity, 

Ep (-k^ F(kt'^)) = M.,iO E^ (-r A;" F (k (i^)")) dC 

Applying the inverse Fourier transformation we get the first fundamental 
solution G^^, 

G^.Jt^x) = M.iOGl^ {t\x/0 J (50) 

or, equivalently 

G';At. ^) = £ M,{0 Gl^ (1, y/0 ^ (51) 

where y = x/f . Similarly the second fundamental solution is given by the 
Mellin convolution 

^) = Gi',„ (1, y/0 J (52) 

6 Concluding remarks. 

Wave equations with fractional-order Laplacians and their non-local generaliza- 
tions allow for a new class of models compatible with very general anisotropy. 
They are compatible with the principles of viscoelasticity. 

Fairly explicit expressions have been obtained for the fundamental functions 
of the Cauchy problem. While solutions by the inverse Laplace and inverse 
Fourier transformation are intractable, solutions in the form of a Mellin convo- 
lution are easier to obtain and analyze. Analysis of the attenuation and dis- 
persion is however more difficult than in local viscoelasticity. For a space-time 
fractional equation the dispersion and attenuation are defined by an equation 
pp'^ + M|k|" = 0. It is easy to see that |k| ~oo Aui'^ for large w := ip. Thus 
the attenuation grows at a sublinear rate in the high-frequency region. 
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A Some definitions 

A function / on the open positive real half- line ]0,cx)[ is said to be completely 
monotone (CM) if it is infinitely differentiable and (— 1)"D" f{t) > for every 
non- negative integer n and t > 0. 

A CM function can have a singularity at 0. A CM function is locally inte- 
grable (LICM) if it its integral over the segment [0, 1] is finite. 

A measurable real function / on [0, cxd[ is said to be causally positive definite 
(CPD) if 

f{y)if{x~y)ip{x) > 







for every compactly s upported real test functi on ip. 

By a theorem in iGripenberg et al" 199d | every LICM function is CPD. 



Gripenberg's extension of Bochner's theorem asserts that every continuous CPD 
function is equal for x > to the Fourier transform of a finite positive Radon 
measure. For our purposes a the restriction to bounded functions is unwelcome. 
The Bochner-Schwartz theorem is an extension of Bochner's theorem to positive 
definite tempered distributions: 

Theorem A.l A causal positive definite tempered distribution f is the Fourier 
transform of a positive tempered Radon measure: 



/oo 
e~*-^/x(dk), a;>0 
-oo 



where is a tempered Radon measure (i.e. a non-negative tempered distribu- 
tion). This theorem captures functions with singularities typical of Gelfand's 
homogeneous distributions. 

B Wright functions 

The Wright function depend on two parameters. They are defined in terms of 
power series expansion 

oo „ 

W^a,m(^):=E^W^Y— ^' A>-1, A^eC 

n=0 ^ ^' 
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Imp 



Figure 5: The Hankel contour T-L. 



They also have an integral representation 



1 

27ti 



H 



A > -1, At e C 



where T-L denotes the Hankel contour encircling the cut along the negative real 
axis, as shown in Fig. [S) 

The Wright functions Af^ and are special cases of the Wright function 



(53) 
(54) 



The best review paper on the Wright function and is iMainardil [201 ij . The 
function N~f has not been studied before. Their series expansions are 



M,{z) = 



i-zY 



-'n!r(l-7(n + l) 



-j^ oo - \n—l 

1 ) 7T (n — 1)! 

-z)" 



7V,(z) = J2 



^-n!r(2-7(n+l) 
These functions also have the integral representations 



1 



27Ti 



(56) 

(57) 
(58) 



Num erical methods of ev aluating the Wright functions can be found in lLuchko et al.l 
2010l| andlLudiko! |2008 |. 
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The Mellin transform of M^{z) is r(l - s)/r(l - 7s) and 

1 e+ioo Yd - s] 

MJz) ^ — / -7^ z-'ds for z > 0, < 7 < 1, < £ < 1 (59) 

2™7-£-ioo r(l-7s) 

As a rough check note that for z > 1 the Bromwieh contour can be closed in 
the right s half-plane. The contribution of the poles of r(l — s) at s = n + 1, 
n = 0, 1, 2, ... is 

' n!r(l-(n+l)7) ' 
Similarly, the Mellin transform of N^{z) is r(l - s)/r(2 - 7s). 



C Proof of identity (IlSj) and related results. 

Lemma C.l 

Proof. Since r(s)r(l — s) — n/ sin(7Ts), the integrand on the right-hand side 
of equation (|60)) has poles at s = n = 0, ±1, ±2, . . .. Stirling's formula 



r(z) = V2nz-'^/^{z/ey{l + 0[l/z]) for |arg(z)| <7t-e, e>0 
implies that for Res — > — 00 

z-'/T{l^l3s) ^,^00 h^sy/yV2^ exp ((/3s - 1) [ln(l - ps) - 1] - s ln(z)) 

The first term in the exponent dominates and therefore the integrand vanishes 
for Res —00 faster than l/|s|. Hence the Bromwieh contour can be closed by 
a half-circle at infinity in the left half of the complex s-plane without changing 
the value of the integral. The only singularities inside the closed contour thus 
obtained are the poles at s = 0,-1,-2,. Their residues are (— z)"/r(l -I- /3n). 
Hence the right-hand side of equation (pO)) is equal to 



-j-zy_ 

^ r(l -hn/3) 



I I I -t- 7/, n 1 

n,=0 



The proof of equation (IBTj) is analogous. □ 
The integrands on the right hand side of ((60|) and ([6T|) are thus Mellin transforms 
of the left-hand sides: 



p ( \ s-x. r(s)r(i-s) 

Efi^2{.-x)x ds = _ (63) 
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Lemma C.2 



M^iy) = 7^, I 7^7^^T^^2/"' els < 7 < 1 (64) 



NJy) = — f —, -v-" As < 7 < 1 (65) 

^^^^ 27ti7gr(2 + 7(s-l))^ ' ^ ' 

Proof. The proof is similar to the proof of Lemma [C.ll The residues of r(s) 
at s = — n, n = 0, 1, . . .. are (— 1)"/*^!: hence the integral on the right-hand side 
of (|64)) is equal to 



^^^n!r(l-7(n+l))' 

which is equal to the function on the left-hand side. 

The other equation is proved in the same way. □ 

We are now ready to prove the following important relations 

/•CO 

Efi i-K") ^ / M^iO E„ (-«)") de (66) 
Jo 

Ep,2 (-/^") = / N^iO Eo. (-«)") d^ (67) 
Jo 

where < 7 = (3/ a < 1. 

Equation (|60p implies that 

'^•<-"°' = 2^/^T^''-^' (68) 

2ma r(l - 7s) 



while equation (j64p implies that 



We now recall that the Mellin transform of the Mellin convolution 

of two functions / and g is equal to the product of the Mellin transforms of / 
and g. Hence 

/•OO -I 1> /"OO 

Ep = / jAMl/O E^ (-(/^/O") 4 = / ^^7(0 Eo. (-«)") dC 

Jo t, Jo 

This proves equation (j66p . Equation (1671) is proved in a similar way. 

The identities (l66l) and (|67)) have an important corollary. Applying the 
inverse Fourier transform to both sides of either identity and noting p6p we 



20 



obtain the relations 



G/3,a(l,a;) 



/■°° df 

/ Af^(OGo,a(l,x/0^ 

Jo s 



iV7(0Ga,o(l,x/0 , 
S 



(70) 
(71) 



D The functions M2/3{z) and A^2/3(2;)- 

A method for num e rical computat ion of the Wright function can be found in 
Luchko et all (2010( . iLuchkol |2008 |. 

Explicit repr esentations of exist for 7 = 1/3, 1/2 and 2/3. The first two 
can be found in iMainardil 201 1| . Since we are interested in wave equations, 
only 7 > 1/2 is of interest to us. 

The function M2/h (z) can be expressed in t e rms o f the Airy function using 
the results obtained in iHanyga and Seredvnskal 2002l | for the function 



/P(i,Ai,A2) 



9 H3) 

d\2 



(i,Ai,A2) 



where 



/i^"* (^: Al, A2) 



1 

27Ti 



(3) 

and B denotes the Bromwich contour. The function /2 can be expressed in 
terms of the Airy function as follows: 



(72) 



(73) 



32/3 



(74) 

On the other hand straightforward transformations of the integration vari- 
able in (|73)) yield the relation 



Al 



2Af 



,/3(z) = /P) (.-^/M,0) /z 



(75) 



It is to be noted that the Hankel contour in (f57| can be replaced by the Bromwich 
contour. 

Equation (|75p can be worked out in terms of the Airy function 



/3(z) = [3-1/3 2 Ai (^73^/3^ - 3i/3Ai' (zV3'/') 



Mo 



or in terms of the modified Bessel function of the second kind 



-22^/27 



-22^/27 



(76) 



(77) 
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A simple expression for the function t/^/"^'") (y) can be derived from equa- 
tion (|77| : 

Ui'/'''\y) = ^ J^^ Z{0X^{2'/'y/ (se^^)) dC (78) 

where 

Z(C) := (Xi/3(C)+i^2/3(C)) c-'^ (79) 

In order to express N2/3 in terms of the Airy functions, we define a new 
function 

F{t, Ai, X,) := f^'\t, e, X,) ^ ^ s-^^' e^ e"^^ ^^^^-^..^3 
In terms of this function 

Substituting (j74l) the following alternative formula is derived 

iV2/3(z) = 32/^z-i£c^/2Ai (3-4/3 C^) cxp(^-Ac') dC (80) 
N2/3 can be further simplified by substituting 

Ai(x) = VV2 ,1/2^^/3 (^^x3/2^ 

which yields the formula 

o /■2z^/27 

^2/3(^) = ^^7^y^ y~'/'i^i/3(2/)e-My (81) 

or 

7V2/3(^) = ^[r(l/3) 2F2 (1/6, 1/2; 1/3, 1/2; -4zV27) 
Ttvo 

- yr(2/3) 2F2 (5/6, 7/6; 5/3, 13/3; -4z727)] (82) 
The functions M2/3 and A^2/3 are shown in Fig.[S] 

E The inverse Fourier transform of E'^^ (— 

Applying the inverse Fourier transformation to both sides of the identity 

^-^ 11 + an 

n— ^ ^ 
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Figure 6: Plots of M2/^{z) (solid line) and ^''2/3(2;) (dashed line). 



and recalling the Fourier t ransform of |fc|" is — 2 sin(na7T/2) r(l + an) \x\ " 
Gel'fand and Shilovl 1964l | we have 



— / e"=^£;„(-|fc|") dfc = -2Ini^(-l)"|x|-""-i (^e'""/2j 



2 T 1 

■; — r-lni ; — ; -. — 77- 

\x\ l+a;-"ei""/2 



so that 

G„,„(l,x)=X„(x) = — / e''=-i?„.„(-|fc|") dfc 

sin(7Ta/2) 



IQ-I 



7T + 2|a;|" cos(7ta/2) + 1 



(83) 



Mainardi et alj |2001 |. 

The second Green's function Ha a is somewhat more difficult to calculate 

Ha,M,x) = Yaix) = — e'^-Ea^^Hkn dfc ^^(|x|) 



But 



Hence 



r(2 + an) 



9 / 
F'(y) = --Im^(-1)" ( 



ji=0 



2sin(a7t/2) 



y 



a-2 



y2a _)_ 2yO' C0s(a7T/2) + 1 
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and F{y) vanishes at infinity. Hence 



/■oo a — 2 

F(y) = -2 sm( an/ 2) — - — -, dz (84) 

V / ; ^2a 2z" cos(a7T/2) + 1 

Unlike Ga,a, the function Ha^a is not an algebraic function. For a = 1/2 it 
can be expressed in terms of the logarithm. 

Combining equations dHU), dTO]), ^ and dSH]) we get fairly exphcit 

integral representations of the fundamental solutions of the Cauchy problem in 
one and three dimensions. 
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